A parallel solver for reaction-diffusion systems in computational electrocardiology

In this work, a parallel three-dimensional solver for numerical simulations in computational electrocardiology is introduced and studied. The solver is based on the anisotropic Bidomain %(AB) cardiac model, consisting of a system of two degenerate parabolic reaction-diffusion equations describing the intra and extracellular potentials of the myocardial tissue. This model includes intramural fiber rotation and anisotropic conductivity coefficients that can be fully orthotropic or axially symmetric around the fiber direction. %In case of equal anisotropy ratio, this system reduces to The solver also includes the simpler anisotropic Monodomain model, consisting of only one reaction-diffusion equation. These cardiac models are coupled with a membrane model for the ionic currents, consisting of a system of ordinary differential equations that can vary from the simple FitzHugh-Nagumo (FHN) model to the more complex phase-I Luo-Rudy model (LR1). The solver employs structured isoparametric $Q_1$ finite elements in space and a semi-implicit adaptive method in time. Parallelization and portability are based on the PETSc parallel library. Large-scale computations with up to $O(10^7)$ unknowns have been run on parallel computers, simulating excitation and repolarization phenomena in three-dimensional domains

Parallel multilevel solvers for the cardiac electro-mechanical coupling

We develop a parallel solver for the cardiac electro-mechanical coupling. The electric model consists of two non-linear parabolic partial differential equations (PDEs), the so-called Bidomain model, which describes the spread of the electric impulse in the heart muscle. The two PDEs are coupled with a non-linear elastic model, where the myocardium is considered as a nearly-incompressible transversely isotropic hyperelastic material. The discretization of the whole electro-mechanical model is performed by Q1 finite elements in space and a semi-implicit finite difference scheme in time. This approximation strategy yields at each time step the solution of a large scale ill-conditioned linear system deriving from the discretization of the Bidomain model and a non-linear system deriving from the discretization of the finite elasticity model. The parallel solver developed consists of solving the linear system with the Conjugate Gradient method, preconditioned by a Multilevel Schwarz preconditioner, and the non-linear system with a Newton\u2013Krylov-Algebraic Multigrid solver. Three-dimensional parallel numerical tests on a Linux cluster show that the parallel solver proposed is scalable and robust with respect to the domain deformations induced by the cardiac contraction

A comparison of coupled and uncoupled solvers for the cardiac Bidomain model

The aim of this work is to compare a new uncoupled solver for the cardiac Bidomain model with a usual coupled solver. The Bidomain model describes the bioelectric activity of the cardiac tissue and consists of a system of a non-linear parabolic reaction-diffusion partial differential equation (PDE) and an elliptic linear PDE. This system models at macroscopic level the evolution of the transmembrane and extracellular electric potentials of the anisotropic cardiac tissue. The evolution equation is coupled through the non-linear reaction term with a stiff system of ordinary differential equations (ODEs), the so-called membrane model, describing the ionic currents through the cellular membrane. A novel uncoupled solver for the Bidomain system is here introduced, based on solving twice the parabolic PDE and once the elliptic PDE at each time step, and it is compared with a usual coupled solver. Three-dimensional numerical tests have been performed in order to show that the proposed uncoupled method has the same accuracy of the coupled strategy. Parallel numerical tests on structured meshes have also shown that the uncoupled technique is as scalable as the coupled one. Moreover, the conjugate gradient method preconditioned by Multilevel Hybrid Schwarz preconditioners converges faster for the linear systems deriving from the uncoupled method than from the coupled one. Finally, in all parallel numerical tests considered, the uncoupled technique proposed is always about two or three times faster than the coupled approach

Role of infarct scar dimensions, border zone repolarization properties and anisotropy in the origin and maintenance of cardiac reentry

Cardiac ventricular tachycardia (VT) is a life-threatening arrhythmia consisting of a well organized structure of reentrant electrical excitation pathways. Understanding the generation and maintenance of the reentrant mechanisms, which lead to the onset of VT induced by premature beats in presence of infarct scar, is one of the most important issues in current electrocardiology. We investigate, by means of numerical simulations, the role of infarct scar dimension, repolarization properties and anisotropic fiber structure of scar tissue border zone (BZ) in the genesis of VT. The simulations are based on the Bidomain model, a reaction-diffusion system of Partial Differential Equations, discretized by finite elements in space and implicit-explicit finite differences in time. The computational domain adopted is an idealized left ventricle affected by an infarct scar extending transmurally. We consider two different scenarios: i) the scar region extends along the entire transmural wall thickness, from endocardium to epicardium, with the exception of a BZ region shaped as a central sub-epicardial channel (CBZ); ii) the scar region extends transmurally along the ventricular wall, from endocardium to a sub-epicardial surface, and is surrounded by a BZ region (EBZ). In CBZ simulations, the results have shown that: i) the scar extent is a crucial element for the genesis of reentry; ii) the repolarization properties of the CBZ, in particular the reduction of IKs and IKr currents, play an important role in the genesis of reentrant VT. In EBZ simulations, since the possible reentrant pathway is not assigned a-priori, we investigate in depth where the entry and exit sites of the cycle of reentry are located and how the functional channel of reentry develops. The results have shown that: i) the interplay between the epicardial anisotropic fiber structure and the EBZ shape strongly affects the propensity that an endocardial premature stimulus generates a cycle of reentry; ii) reentrant pathways always develop along the epicardial fiber direction; iii) very thin EBZs rather than thick EBZs facilitate the onset of cycles of reentry; iv) the sustainability of cycles of reentry depends on the endocardial stimulation site and on the interplay between the epicardial breakthrough site, local fiber direction and BZ rim

A clinical-in silico study on the effectiveness of multipoint bicathodic and cathodic-anodal pacing in cardiac resynchronization therapy.

Up to one-third of patients undergoing cardiac resynchronization therapy (CRT) are nonresponders. Multipoint bicathodic and cathodic-anodal left ventricle (LV) stimulations could overcome this clinical challenge, but their effectiveness remains controversial. Here we evaluate the performance of such stimulations through both in vivo and in silico experiments, the latter based on computer electromechanical modeling. Seven patients, all candidates for CRT, received a quadripolar LV lead. Four stimulations were tested: right ventricular (RVS); conventional single point biventricular (S-BS); multipoint biventricular bicathodic (CC-BS) and multipoint biventricular cathodic-anodal (CA-BS). The following parameters were processed: QRS duration; maximal time derivative of arterial pressure (dPdtmax); systolic arterial pressure (Psys); and stroke volume (SV). Echocardiographic data of each patient were then obtained to create an LV geometric model. Numerical simulations were based on a strongly coupled Bidomain electromechanical coupling model. Considering the in vivo parameters, when comparing S-BS to RVS, there was no significant decrease in SV (from 45 ± 11 to 44 ± 20 ml) and 6% and 4% increases of dPdtmax and Psys, respectively. Focusing on in silico parameters, with respect to RVS, S-BS exhibited a significant increase of SV, dPdtmax and Psys. Neither the in vivo nor in silico results showed any significant hemodynamic and electrical difference among S-BS, CC-BS and CA-BS configurations. These results show that CC-BS and CA-BS yield a comparable CRT performance, but they do not always yield improvement in terms of hemodynamic parameters with respect to S-BS. The computational results confirmed the in vivo observations, thus providing theoretical support to the clinical experiments

Incorporating Inductances in Tissue-Scale Models of Cardiac Electrophysiology

In standard models of cardiac electrophysiology, including the bidomain and monodomain models, local perturbations can propagate at infinite speed. We address this unrealistic property by developing a hyperbolic bidomain model that is based on a generalization of Ohm's law with a Cattaneo-type model for the fluxes. Further, we obtain a hyperbolic monodomain model in the case that the intracellular and extracellular conductivity tensors have the same anisotropy ratio. In one spatial dimension, the hyperbolic monodomain model is equivalent to a cable model that includes axial inductances, and the relaxation times of the Cattaneo fluxes are strictly related to these inductances. A purely linear analysis shows that the inductances are negligible, but models of cardiac electrophysiology are highly nonlinear, and linear predictions may not capture the fully nonlinear dynamics. In fact, contrary to the linear analysis, we show that for simple nonlinear ionic models, an increase in conduction velocity is obtained for small and moderate values of the relaxation time. A similar behavior is also demonstrated with biophysically detailed ionic models. Using the Fenton-Karma model along with a low-order finite element spatial discretization, we numerically analyze differences between the standard monodomain model and the hyperbolic monodomain model. In a simple benchmark test, we show that the propagation of the action potential is strongly influenced by the alignment of the fibers with respect to the mesh in both the parabolic and hyperbolic models when using relatively coarse spatial discretizations. Accurate predictions of the conduction velocity require computational mesh spacings on the order of a single cardiac cell. We also compare the two formulations in the case of spiral break up and atrial fibrillation in an anatomically detailed model of the left atrium, and [...].Comment: 20 pages, 12 figure

Cardiac anisotropy in boundary-element models for the electrocardiogram

The boundary-element method (BEM) is widely used for electrocardiogram (ECG) simulation. Its major disadvantage is its perceived inability to deal with the anisotropic electric conductivity of the myocardial interstitium, which led researchers to represent only intracellular anisotropy or neglect anisotropy altogether. We computed ECGs with a BEM model based on dipole sources that accounted for a “compound” anisotropy ratio. The ECGs were compared with those computed by a finite-difference model, in which intracellular and interstitial anisotropy could be represented without compromise. For a given set of conductivities, we always found a compound anisotropy value that led to acceptable differences between BEM and finite-difference results. In contrast, a fully isotropic model produced unacceptably large differences. A model that accounted only for intracellular anisotropy showed intermediate performance. We conclude that using a compound anisotropy ratio allows BEM-based ECG models to more accurately represent both anisotropies

Optimized schwarz methods and model adaptivity in electrocardiology simulations

[No abstract available

Comparison of diffusion tensor imaging by cardiovascular magnetic resonance and gadolinium enhanced 3D image intensity approaches to investigation of structural anisotropy in explanted rat hearts

Background: Cardiovascular magnetic resonance (CMR) can through the two methods 3D FLASH and diffusion tensor imaging (DTI) give complementary information on the local orientations of cardiomyocytes and their laminar arrays. Methods: Eight explanted rat hearts were perfused with Gd-DTPA contrast agent and fixative and imaged in a 9.4T magnet by two types of acquisition: 3D fast low angle shot (FLASH) imaging, voxels 50 × 50 × 50 μm, and 3D spin echo DTI with monopolar diffusion gradients of 3.6 ms duration at 11.5 ms separation, voxels 200 × 200 × 200 μm. The sensitivity of each approach to imaging parameters was explored. Results:The FLASH data showed laminar alignments of voxels with high signal, in keeping with the presumed predominance of contrast in the interstices between sheetlets. It was analysed, using structure-tensor (ST) analysis, to determine the most (v 1 ST ), intermediate (v 2 ST ) and least (v 3 ST ) extended orthogonal directions of signal continuity. The DTI data was analysed to determine the most (e 1 DTI ), intermediate (e 2 DTI ) and least (e 3 DTI ) orthogonal eigenvectors of extent of diffusion. The correspondence between the FLASH and DTI methods was measured and appraised. The most extended direction of FLASH signal (v 1 ST ) agreed well with that of diffusion (e 1 DTI ) throughout the left ventricle (representative discrepancy in the septum of 13.3 ± 6.7°: median ± absolute deviation) and both were in keeping with the expected local orientations of the long-axis of cardiomyocytes. However, the orientation of the least directions of FLASH signal continuity (v 3 ST ) and diffusion (e 3 ST ) showed greater discrepancies of up to 27.9 ± 17.4°. Both FLASH (v 3 ST ) and DTI (e 3 DTI ) where compared to directly measured laminar arrays in the FLASH images. For FLASH the discrepancy between the structure-tensor calculated v 3 ST and the directly measured FLASH laminar array normal was of 9 ± 7° for the lateral wall and 7 ± 9° for the septum (median ± inter quartile range), and for DTI the discrepancy between the calculated v 3 DTI and the directly measured FLASH laminar array normal was 22 ± 14° and 61 ± 53.4°. DTI was relatively insensitive to the number of diffusion directions and to time up to 72 hours post fixation, but was moderately affected by b-value (which was scaled by modifying diffusion gradient pulse strength with fixed gradient pulse separation). Optimal DTI parameters were b = 1000 mm/s2 and 12 diffusion directions. FLASH acquisitions were relatively insensitive to the image processing parameters explored. Conclusions: We show that ST analysis of FLASH is a useful and accurate tool in the measurement of cardiac microstructure. While both FLASH and the DTI approaches appear promising for mapping of the alignments of myocytes throughout myocardium, marked discrepancies between the cross myocyte anisotropies deduced from each method call for consideration of their respective limitations